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Abstract 

In  order  to  formulate  mathematical  conjectures  likely  to  be  true,  a 
number  of  base  cases  must  be  determined.  However,  many  combinato¬ 
rial  problems  are  NP-hard  and  the  computational  complexity  makes  this 
research  approach  difficult  using  a  standard  brute  force  approach  on  a  typ¬ 
ical  computer.  One  sample  problem  explored  is  that  of  finding  a  minimum 
identifying  code.  To  work  around  the  computational  issues,  a  variety  of 
methods  are  explored  and  consist  of  a  parallel  computing  approach  using 
Matlab,  a  quantum  annealing  approach  using  the  D-Wave  computer,  and 
lastly  using  satisfiability  modulo  theory  (SMT)  and  corresponding  SMT 
solvers.  Each  of  these  methods  requires  the  problem  to  be  formulated  in 
a  unique  manner.  In  this  paper,  we  address  the  challenges  of  computing 
solutions  to  this  NP-hard  problem  with  respect  to  each  of  these  methods. 


1  Problem  Statement  and  Background 

First  introduced  in  1998  [4],  an  identifying  code  for  a  graph  G  is  a  subset  of 
the  vertices,  S  C  V(G),  such  that  for  each  v  G  V{G)  the  subset  of  vertices 
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Figure  1:  The  2-ary  de  Bruijn  graph  of  order  3,  or  13(2,  3). 

of  S  that  are  adjacent  to  v  is  non-empty  and  unique.  That  is,  each  vertex  of 
the  graph  is  uniquely  identifiable  by  the  non-empty  subset  of  vertices  of  S  to 
which  it  is  adjacent.  More  formally,  let  N(v)  be  the  set  of  vertices  adjacent  to 
v  and  B(v )  =  N(v)  U  {u}.  Then  we  require  that  any  two  vertices  u,v  £  V(G) 
have  different  identifying  sets,  or  more  precisely  that  we  must  have  SnB(v)  ^ 
SnB(u),  and  also  that  both  SDB(v),  SC\B(u)  ^  0.  The  combinatorial  problem 
of  finding  minimum  identifying  codes  has  been  shown  to  be  NP-complete  [3], 
but  also  has  many  potential  real-world  applications. 

As  some  NP-complete  problems  are  notorious  for  having  special  graph  classes 
on  which  there  are  simple  solutions,  previous  research  has  focused  on  the  class 
of  de  Bruijn  networks  [2].  This  paper  explores  the  problem  of  finding  the  mini¬ 
mum  size  of  an  identifying  code  over  the  undirected  de  Bruijn  graph  using  three 
different  methods.  Section [2]  explores  a  parallel  computing  algorithm,  Section  [3] 
discusses  mapping  the  problem  onto  the  D-Wave  computer  and  using  a  quan¬ 
tum  annealing  approach,  and  finally  Section  [f]  illustrates  the  method  of  using 
Satisfiability  Modulo  Theory  (SMT)  solvers. 

While  many  of  our  examples  and  data  revolve  around  the  class  of  de  Bruijn 
graphs,  the  methods  discussed  throughout  can  easily  be  applied  to  arbitrary 
graphs.  For  reference,  we  provide  some  of  the  basic  definitions  and  background 
here.  The  undirected  d-ary  de  Bruijn  graph  of  order  n,  denoted  B(d,n),  is  the 
graph  with  the  following  vertex  and  edge  sets. 

V  =  {x^.-.Xn  I  Xi  e  {0, 1,  1}} 

E  =  {(x,y)  \  x2x3  . .  .xn  =  yiy2  ■  ■  .yn-i  ov  xix2  ■  ■ -xn-i  =  y2y3  ■  ■ -yn} 

For  example,  the  graph  B( 2,  3)  is  illustrated  in  Figure  [l]  These  graphs  have 
many  useful  properties  for  applications,  such  as  having  a  relatively  high  number 
of  nodes,  a  low  degree  at  each  node,  and  many  short  paths  between  any  two 
nodes. 


2  Parallel  Computing 

The  Information  Directorate  of  the  Air  Force  Research  Laboratory  is  home 
to  one  of  the  Department  of  Defense’s  high-performance  computing  systems 
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(HPC).  Installed  on  the  HPC  is  Matlab,  including  the  Parallel  Computing  Tool¬ 
box.  The  pseudocode  in  Algorithm  [l]  describes  our  brute  force  algorithm,  im¬ 
plemented  in  Matlab. 


Algorithm  1  Brute  Force  Algorithm 

1:  procedure  BruteForce((I,  n,  k)  >  B(d,n)  and  subset  size  k 

2:  Create  list  of  all  subsets  of  k  nodes 

3:  for  i  =  1  :  (dfc )  do 

4:  if  subset  i  is  a  valid  identifying  code  then 

5:  Display  subset  i  to  user 

6:  else 

7:  Do  nothing 

8:  end  if 

9:  end  for 

10:  end  procedure 


Parallelizing  our  algorithm  takes  two  steps.  The  first  step  is  to  replace  “For” 
on  line  3  with  “Parfor” .  This  indicates  to  Matlab  to  use  the  parallel  computing 
toolbox  and  run  each  loop  iteration  independently.  The  second  step  requires 
moving  the  construction  of  subsets  inside  the  Parfor  loop.  Because  of  the  expo¬ 
nential  increase  in  the  number  of  subsets  created,  it  is  more  efficient  to  generate 
each  subset  within  the  loop  and  discard  it  after  the  iteration  than  to  store  all 
(' dk  )  fc-subsets  and  traverse  through  the  list.  This  is  done  using  a  /e-subset  un¬ 
ranking  algorithm.  Two  of  these  algorithms  (from  [5])  are  listed  as  Algorithms 
[2]  and  [3]  These  unranking  functions  allow  us  to  completely  parallelize  the  brute 
force  algorithm,  and  the  results  obtained  are  listed  in  Figure  [2] 


Algorithm  2  Revolving  Door  Unranking  Algorithm 


procedure  REVDoOR(r,  k,  n) 
x  =  n 

for  i  =  k  :  1  do 

while  >  r  do 

x  =  x  —  1 

end  while 

ti  =  x  +  1 


>  subset  index,  subset  size,  set  size 


8 

9 

10 

11 


end  for 

end  procedure 
return  T  =  (fi,f2, . . .  ,4) 


3  D-Wave  Quantum  Annealing  Machine 

Under  the  collaborative  effort  “Adiabatic  Quantum  Computing  Applications 
Research”  (14-RI-CRADA-02)  between  the  Information  Directorate  and  Lock- 
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Algorithm  3  Lexicographic  Unranking  Algorithm 


1 

procedure  LexUnrank/t,  k,  n)  >  subset  index,  subset  size,  set  size 

2 

X  =  1 

3 

for  i  =  1  :  k  do 

4 

while  r  >  (nkZ*)  do 

5 

II 

1 

S 

1  1 
*»•  H 

6 

x  =  x  +  1 

7 

end  while 

8 

ti  =  X 

9 

X  =  X  +  1 

10 

end  for 

11 

end  procedure 

12 

return  T  =  (fi,f2, .  -  -  ,tk) 

d\n 

2  3  4  5 

2 

x  4  6  12 

3 

4  9 

4 

5 

5 
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Figure  2:  Results  for  B(d,n)  obtained  using  HPC 

heed  Martin  Corporation,  we  aim  to  extend  the  results  obtained  by  the  HPC 
system.  In  general,  the  D-Wave  machine  can  address  a  class  of  Ising  problems 
natively  by  the  hardware.  As  stated  in  the  D-Wave  user  documents,  “The  D- 
Wave  hardware  can  be  viewed  as  a  hardware  heuristic  which  minimizes  Ising 
objective  functions  using  a  physically  realized  version  of  quantum  annealing.” 
[7]  The  Ising  model  is  an  energy  minimization  problem  of  -1/+1- valued  vari¬ 
ables.  It  can  be  converted  to  a  quadratic  unconstrained  binary  optimization 
(QUBO)  problem  that  uses  0/1-valued  variables,  and  so  they  are  often  used 
interchangeably. 

3.1  Binary  Optimization  Model 

We  present  a  binary  optimization  formula  for  the  minimum  identifying  code 
problem.  Adjustments  must  be  made  to  create  a  quadratic  version.  We  will 
define  this  model  using  three  separate  functions:  one  to  show  that  the  set  has 
the  correct  size,  one  to  show  that  the  set  is  dominating,  and  one  to  show  that 
the  set  is  separating  (or  identifying). 

3.1.1  Variable  Definitions 

We  will  use  the  notation  B{y)  for  v  £  V(G),  where  B(v)  =  N(v)U{v}.  In  other 
words,  B(y)  is  the  set  containing  all  vertices  adjacent  to  v,  plus  v  itself.  This  is 
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referred  to  in  graph  theory  as  the  ball  of  radius  one  centered  at  v. 

We  define  the  variables  as  follows. 

=  /  !-  ^  i  e  B(v); 

Vl  \  0,  otherwise. 

3.1.2  Set  S  has  size  k 

We  define  the  first  function,  Ha ,  as  follows. 

Ha  [k  xw) 

0  iff  |S|  =  k. 

3.1.3  Set  S'  is  a  dominating  set 

By  definition,  this  is  equal  to  \/v  G  G,  B(v)  D  S  ^  0.  This  is  equivalent  to  the 
following. 


\/v  G  G,  B(v)  n  S  7^  0  o  ( xvv  =  1)  V  (  ^  xuv  >  1  j 

\uv£E  / 

(1  Xvv  =  0)  V  ~ 1  (  ^  ^  XUv  =  0  j 
\uv£E  / 

o  (i-iOT=o)v  n  (i ^ *««) = ° ) 

\uvEE  / 

From  this  statement,  we  get  the  following  equation  for  our  second  function. 

Hb  —  —  xvv)  '  —  xuv)) 


3.1.4  Set  S  is  a  separating  set 

By  definition,  this  is  equal  to  Vx,y  G  G,  ( B{x )  D  S)A(B(y)  fl  S)  7^  0.  This  is 
equivalent  to  the  following  for  a  specific  pair  x  ^  y. 


(B{x)  (~l  S)A(B(y)  n  S)  ^  0  o  G  {B(x)  n  S)A(B(y)  n  S) 

G->-  3v,  ( v  G  B{x)  fl  S)  ©  (w  G  B{y)  fl  S) 

G7  dl?,  ( Xxv  1)  ©  ( Xyy  —  1) 

G7  dl?,  (1  ( Xxv  ©  Xyy )  0) 

G7  |  |  (1  Xxv  Xyy  )  0 

V 

From  this  statement,  we  get  the  following  equation  for  our  third  function, 
summed  over  all  pairs  x,  y. 

He  =  Ylx  'Yhy^x  rii>(l  ~  Xxv  —  Xyy) 
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3.1.5  The  Binary  Optimization  Model 

From  these  three  functions,  our  binary  optimization  model  is  the  following. 

H(S)  =  Ha(S)  +  Hb(S)  +  HC(S) 

=  0  iff  S  is  an  identifying  code. 

Note  that  while  this  does  provide  a  binary  optimization  model  for  our  prob¬ 
lem,  it  is  not  quadratic.  In  order  to  convert  H(S)  to  a  quadratic  binary  equation, 
each  higher  order  term  must  be  replaced  with  several  new  variables.  While  this 
is  possible,  it  is  a  time-consuming  and  arduous  process  that  introduces  many 
new  variables.  Hence  this  approach  will  likely  not  be  the  most  efficient  imple¬ 
mentation. 


3.2  Integer  Program  Formulation 

From  [3],  we  have  the  following  integer  program  formulation  of  the  minimum 
identifying  code  problem  in  graphs. 

First,  we  define  the  modified  adjacency  matrix  as  follows.  It  is  the  adjacency 
matrix  plus  the  identity  matrix. 

.  f  t,  if  (i,j)  €  E  or  i  =  j; 

]^0,  otherwise. 

Using  this  definition,  we  see  that  a  ball  of  radius  1  centered  at  vertex  i  is 
given  by  the  following  vector. 

B(l)  =  A‘2j ,  .  .  ■  , 


Our  vertex  subset  S  is  defined  as  the  following  vector. 


S  =  [si,  s2,  •  •  • ,  sn]T  where  S; 


1,  if  i  £  S; 

0,  otherwise. 


To  compare  two  identifying  sets  with  respect  to  S  for  vertices  i  and  j .  the 
following  expression  computes  the  size  of  (B(i)  fl  S)A(B(j)  n  S). 


n 

^  ^  \A-ki  -^kj  |  '  Sk 
fc= 1 


This  implies  that  in  order  for  S  to  be  a  valid  identifying  code,  we  must  have 
the  following  inequality  satisfied  for  all  pairs  of  vertices  i  and  j. 

n 

~  |  Aki  Akj  |  *  Sk  ^  1 
fc= 1 


For  the  dominating  property  to  be  satisfied,  we  require  the  following  addi¬ 
tional  inequality. 


A-S>  1T 
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Thus  our  integer  program  is  given  by  the  following, 
min  |  S' | 

s-t.  Ylk= 1 1 Aki  -  Akj\  ■  sk  >  1,  Vi  7^  j 

A-S  >  1T 

sk  G  {0, 1} 

In  order  to  use  these  ideas  for  the  D-Wave  machine,  our  constraints  must  be 
equalities.  This  means  we  must  add  binary  slack  variables  for  each  inequality. 
For  the  first  set  of  inequalities,  we  must  determine  an  upper  bound  for  each 
inequality.  Since  these  correspond  to  the  constraint  |(B(i)nS)A(I3(j)nS)|  >  1, 
an  easy  upper  bound  is  given  by  the  following. 

\B(i) |  +  \B(j)\  >  \(B(i )  n  n  S)|  >  1 

For  the  class  of  de  Bruijn  graphs,  we  are  able  to  use  this  to  get  a  bound  on  the 
number  of  slack  variables  needed.  Since  the  maximum  size  of  any  ball  in  B(d,  n) 
is  2d  +  1,  this  gives  us  an  upper  bound  of  size  Ad  +  2  for  this  class  of  graphs. 
Hence  for  each  inequality  in  this  set,  we  must  add  Ad +  2  binary  slack  variables 
to  convert  the  inequality  to  an  equality.  Since  there  are  d  ^ — —  possible  pairs 
i,j,  this  implies  that  we  must  add  a  huge  number  of  binary  slack  variables,  equal 
to  the  following  expression,  just  to  satisfy  the  first  set  of  inequalities. 

dn(dn  —  1)  (2 d  +  1)  slack  variables 

Hence,  this  method  is  not  going  to  be  an  efficient  way  to  map  our  problem  onto 
the  D-Wave  machine. 

3.3  Satisfiability  Formulation 

This  approach  formulates  the  identifying  code  problem  as  a  boolean  satisfiability 
(SAT)  problem.  Each  term  in  the  SAT  problem  is  mapped  to  an  Ising  model. 
The  mapping  introduces  auxiliary  variables  so  that  the  Ising  model  contains 
only  quadratic  and  linear  terms.  A  graph  embedding  technique  is  then  used  to 
map  this  Ising  model  onto  the  connectivity  of  the  D-Wave  chip.  Finally,  gauge 
transformations  are  used  to  mitigate  the  effects  of  intrinsic  control  errors.  We 
will  illustrate  each  step  with  an  example  of  B(d,n)  when  d  =  2  and  n  =  3.  As 
stated  previously,  these  methods  easily  apply  to  any  graph. 

3.3.1  SAT  Formulation 

First,  we  label  the  nodes  of  15(2,3)  from  0  to  dn  —  1  =  7.  Then  we  define  the 
set  of  variables  {x.J  for  i  =  0, 1, . . . ,  7  as  follows. 

_  J  1,  if  node  i  is  included  in  the  identifying  code; 

Xl  (  0,  otherwise. 

Next,  we  look  at  the  ball  for  each  node  and  form  clauses  corresponding  to 
their  domination  constraints. 
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Contents 


Ball 

s(oy 

5(1) 
5(  2) 
B{  3) 
5(4) 
5(5) 
5(6) 
5(7) 


£(000) 

£(001) 

£(010) 

£(011) 

£(100) 

£(101) 

£(110) 

5(111) 


(000,001, 100}  =  (0, 1,4} 

(000,  001,  010,  Oil,  100}  =  {0,1,  2, 3, 4} 
(001,  010, 100, 101}  =  (1,  2,  4,  5} 
(001,011, 101, 110,  111}  =  (1,3,  5, 6,  7} 
(000,001,010, 100,110}  =  (0, 1,2, 4,6} 
(010,  Oil,  101, 110}  =  (2,  3,  5,  6} 

(Oil,  100, 101, 110,  111}  =  (3, 4,  5,  6,  7} 
(Oil,  110, 111}  =  {3,6,7} 


Constraints 

10  V  II  V  X4 

IQ  V  II  V  12  V  13  V  V4 
XI  V  X2  V  X4  V  x$ 

11  V  13  V  15  V  16  v  X 7 
10  V  11  V  12  V  14  V  X6 

12  V  13  V  15  V  XQ 

13  V  14  V  15  V  16  V  Xf 
13  V  16  V  X7 


From  this  set  of  constraints,  we  form  clauses  for  each  pairwise  XOR  of  balls. 
This  is  shown  in  Figure  [3] 

Now  we  can  eliminate  more  specific  clauses  that  are  implied  by  more  general 
clauses.  For  example,  Figure  [4]  shows  which  two-term  constraints  imply  the 
corresponding  larger  constraints.  Hence,  the  only  constraints  that  we  have  left 
are  given  below. 


x2  V  x3 
Xq  V  X2  V  x3 
X2  V  x6 
x0  V  a; 3  V  a;5 
x3  V  x6 
xq  V  £5  V  xq 
X\  V  X2  V  X7 
Xi  V  X4 
X\  V  x3 
X2  V  X4  V  X7 
X2  V  X5  V  X7 
X4  V  X5 

Satisfying  this  set  of  constraints  are  the  four  possible  minimum  solutions, 
given  below. 


{X\,X2,X3,X3} 

{X1,X2,X3,X6\ 

{x2,X3,X4,X3} 

{x2,X4,X3,Xq} 


3.3.2  Mapping  SAT  Clauses  to  Ising  Models 

We  construct  a  Hamiltonian 


U  =  G  Aj})  +  Xj. 

3  i 

Each  of  the  terms  has  the  property  that 


(\J  Xi  is  true 

ieAj 


x*  =  SLTgminHj ({xi,  i  £  A j})  iff 
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We  will  show  momentarily  how  the  T-Lj  are  constructed.  The  last  term  A  >  0  is  a 
penalty  term  that  rewards  smaller  size  codes.  Therefore,  the  minimum  solutions 
(or  ground  states)  of  T~L  are  the  minimum  identifying  codes. 

In  order  to  solve  the  Hamiltonian  using  adiabatic  quantum  optimization, 
we  have  the  further  constraint  that  the  T-L.j  must  contain  only  quadratic  and 
linear  terms  in  the  binary  variables  {aq}.  In  general  to  accomplish  this,  we 
must  introduce  auxiliary  variables,  which  we  will  denote  by  {zi}.  Also,  we  will 
switch  to  the  Ising  model  convention  where  each  of  the  aq  and  Zi  can  take  values 
{+1,  —1}  instead  of  {0, 1}. 

The  mapping  from  OR-clauses  to  Ising  models  that  we  will  use,  namely 
V  Xi  Uj({xui  e  A}), 

i&Aj 

depends  only  on  the  number  of  variables  k  =  \Aj\  in  the  OR-clause.  These 
mappings  for  k  =  2  through  k  =  6  are  represented  diagrammatically  in  Figure 

m 

In  the  diagrams,  numbers  attached  to  a  node  represent  the  linear  coefficients 
in  the  Ising  model,  while  numbers  attached  to  an  edge  represent  the  quadratic 
(coupling)  coefficients  in  the  Ising  model.  For  example,  the  diagram  for  k  =  3 
represents  the  following  Ising  model. 


H3(xi,X2,X3,Zl)  =  XlX2  ~  2xiZl  -  2x2Z\  -  2X2Z\  +  Zix3  +  XX  +  X2  —  3zi  —  x3 

It  can  be  confirmed  that  every  ground  state  of  "^(aq,  X2,  £3,  Zi)  satisfies 
Xi  V  X2  V  X3,  and  conversely  every  combination  of  {aq,  aq,  £3}  that  satisfies 
ii  V  12  V  £3  corresponds  to  a  ground  state  of  7-^3 (aci ,  aq,  aq,  z\). 

3.3.3  Mapping  the  Ising  model  onto  the  D-Wave  processor 

In  general,  the  graph  of  the  Hamiltonian  we  get  from  the  SAT-to-Ising  mapping 
will  not  fit  onto  the  D-Wave  hardware  graph.  The  D-Wave  hardware  graph, 
which  is  called  the  “Chimera”  graph,  is  built  up  of  unit  cells,  each  of  which  is  a 
four  by  four  bipartite  graph,  K dq.  Even  the  simple  Ising  model  for  3-OR  shown 
in  Figure  [5]  cannot  be  mapped  directly  onto  the  D-Wave  hardware  graph.  This 
can  be  seen  from  the  fact  that  our  graph  B( 2,3)  contains  a  3-cycle,  whereas  the 
smallest  cycle  possible  on  the  D-Wave  hardware  graph  is  a  4-cycle. 

Embedding 


Our  first  step  to  embedding  is  to  determine  how  to  map  our  OR-clauses  to 
the  physical  qubits.  One  way  to  embed  the  3-OR  graph  onto  the  D-Wave  is 
shown  in  Figure  [6j  We  have  mapped  the  logical  qubit  z\  to  two  physical  qubits, 
which  are  ferromagnetically  coupled  with  a  coupling  strength  -JFm. 

Similarly,  once  we  construct  the  full  Ising  Hamiltonian  for  the  minimum 
identifying  code  problem,  we  can  use  embedding  to  map  the  graph  of  the  Hamil¬ 
tonian  onto  the  D-Wave  hardware  graph. 
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Figure  6:  Embedding  3-OR  onto  the  Chimera  Graph 


Problem  Decomposition 


If  the  graph  of  the  Hamiltonian  is  too  large  to  embed  onto  our  current  504- 
qubit  D-Wave  hardware  graph,  one  trick  we  can  try  is  to  decompose  the  SAT 
problem  into  smaller  pieces  that  can  be  embedded.  For  example,  in  the  13(2,3) 
example  from  earlier,  one  of  the  terms  is  In  order  for  this  to  be  true,  at 

least  one  of  X2  and  x3  must  be  true.  We  consider  and  solve  each  case  separately. 

If  X2  is  true,  then  so  is  any  OR-clause  containing  X2,  so  we  can  eliminate 
those  from  the  conjunction.  We  are  left  with  the  following  clauses. 

x3  V  x6 
Xo  V  X5  V  Xq 
x\  V  X4 
x\  V  X5 
X4  V  X5 

Similarly,  in  the  case  where  x3  is  true  we  can  eliminate  terms  containing  x3, 
leaving  the  following  clauses. 

x2  V  x6 
Xq  V  x3  V  Xq 
X\  V  X2  V  X7 
X\  V  X4 
X\  V  X5 
X2  V  X4  V  X7 
12  V  15  V  X7 
X4  V  x3 


Both  of  these  subproblems  are  simpler  than  the  original  problem  and  hence 
easier  to  embed.  Whichever  subproblem  yields  the  smaller  identifying  codes  will 
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Figure  7:  Example  Gauge  Transformation 


be  the  solution  of  our  original  problem,  or  if  both  subproblems  have  minimum 
solutions  of  the  same  size,  then  we  can  take  the  union  of  the  two  solution  sets. 

Gauge  Transformations 


Embedding  complex  graphs  leads  to  long  chains ,  i.e.  multiple  physical  qubits 
corresponding  to  the  same  logical  qubit.  In  the  current  D-Wave  embedding 
solver,  all  of  the  physical  qubits  in  a  chain  will  be  ferromagnetically  coupled 
with  some  coupling  strength  -JFm,  which  is  iteratively  increased  to  a  large 
enough  magnitude  to  ensure  that  all  of  the  physical  qubits  in  the  chain  agree 
(most  of  the  time)  and  can  be  treated  as  a  single  logical  qubit. 

However,  due  to  the  characteristics  of  the  D-Wave  design,  embeddings  with 
many  long  ferromagnetically  coupled  chains  will  have  a  greater  tendency  for 
control  precision  errors  which  may  affect  solution  quality.  To  combat  these 
effects,  we  can  utilize  a  gauge  transformation ,  where  we  redefine  a  subset  of  the 
spin  variables  to  be  the  opposite  sign.  Flipping  a  subset  of  the  spins  in  this 
way  induces  a  transformation  on  the  Ising  coefficients:  If  S±  has  been  flipped 
(S[  =  — Si ),  then  h\  =  —hi.  If  one  of  Si  and  S2  has  been  flipped,  then 
J'l2  =  —  Ji2-  But  if  both  of  Si  and  S2  have  been  flipped,  then  J'12  =  Jii- 

Consider  for  example  the  gauge  transformation  shown  in  Figure  [7]  In  the 
figure,  the  red  qubits  are  flipped  by  the  gauge  transformation  while  the  blue 
qubits  are  unchanged.  In  the  first  unit  cell,  all  of  the  horizontal  qubits  are  flipped 
while  the  vertical  qubits  are  unchanged.  In  the  next  unit  cell,  all  of  the  vertical 
qubits  are  flipped  while  the  horizontal  qubits  are  unchanged.  So,  suppose  our 
embedding  contains  the  chain  1  (blue  vertical  qubit),  5  (red  horizontal  qubit),  13 
(blue  horizontal  qubit).  Note  that  each  consecutive  pair  in  the  chain  contains 
exactly  one  flipped  qubit,  so  all  of  the  ferromagnetic  couplings  -JFm  in  the 
chain  will  be  replaced  with  antiferromagnetic  couplings  +JFm  by  the  gauge 
transformation.  By  using  gauge  transformations  like  this,  we  may  be  able  to 
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reduce  the  control  precision  errors  caused  by  embeddings  with  long  chains. 

3.3.4  Results  for  13(2, 4) 

By  combining  all  of  the  above  tricks,  we  obtained  results  for  the  minimum 
identifying  code  problem  on  the  d  =  2,  n  =  4  undirected  de  Bruijn  graph  on  the 
D-Wave  hardware. 

SAT  formulation 


Figure [8] shows  the  full  SAT  formulation  of  the  identifying  code  for  6(2,4). 
This  was  obtained  by  following  the  procedure  outlined  in  section  3.3.1|  It  con¬ 
tains  50  clauses  over  16  variables. 


Problem  Decomposition 


For  the  full  SAT  formulation  in  Figure  [HJ  the  Ising  model  was  too  large  to 
embed  on  the  current  D-Wave  hardware  graph.  The  problem  must  be  decom¬ 
posed  further.  Note  that  two  of  the  terms  are  X3  V  Xu  and  X4  V  X12 .  Thus,  at 
least  one  of  X3  and  X\\  must  be  true,  and  at  least  one  of  X4  and  X12  must  be  true. 
Considering  the  branch  where  X3  and  X4  are  true,  the  SAT  problem  reduces  to 
the  problem  shown  in  Figure  [9]  This  decomposed  formulation  consists  of  just 
24  clauses  over  14  variables. 

SAT  clause  to  Ising  model  mapping 


Using  the  Ising  model  mappings  shown  in  Figure  [5j  we  generated  an  Ising 
model  with  49  auxiliary  variables  {z,;}  for  a  total  of  63  variables.  We  furthermore 
added  the  penalty  term  A  JA  Xi  so  that  the  ground  state  will  be  a  minimum 
identifying  code.  Note  that  this  is  far  better  than  the  1200+  auxiliary  (slack) 
variables  required  in  the  integer  program  formulation  for  this  case. 

Figure  pH]  shows  the  logical  graph  of  the  Ising  model.  In  the  figure,  nodes 
corresponding  to  the  original  14  boolean  variables  are  shown  in  green;  the  re¬ 
maining  nodes  represent  the  auxiliary  variables  added  during  the  SAT-to-Ising 
mapping  process. 

Embedding 


We  used  the  D-Wave  software  function  sapiFindEmbeddingO ,  to  find  an 
embedding  of  this  Ising  model  onto  the  hardware  graph  of  the  504-qubit  D- 
Wave  machine  at  the  USC-Lockheed  Martin  Quantum  Computing  Center  in 
Marina  del  Rey,  CA.  Figure  0  shows  the  embedding  that  we  used,  consisting 
of  253  physical  qubits,  with  a  maximum  chain  length  of  8. 

In  Figure  0  physical  qubits  corresponding  to  the  same  logical  qubit  have 
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Figure  8:  SAT  Formulation  for  6(2,4) 
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Figure  9:  Decomposed  SAT  formulation  for  £3,0:4  true 


Figure  10:  Logical  graph  of  the  Ising  model 
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Figure  11:  Embedding  the  problem  onto  the  D-Wave  hardware 
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the  same  color  and  are  labeled  with  the  same  number.  The  unlabeled  red  qubits 
are  known  faulty  qubits  and  are  not  used. 

Gauge  Transformations 


Since  we  have  an  exact  solution  from  the  HPC  system  on  the  graph  B{ 2, 4), 
we  know  that  the  £3  =  £4  =  1  branch  of  the  problem  should  have  a  ground 
state  with  the  remaining  variables 


Xg  Xiq  X\g  £14  1 

corresponding  to  the  minimum  code  size  of  6.  However,  using  the  standard 
D-Wave  embedding  solver  (which  does  not  yet  support  gauge  transformations 
-  it  uses  all  ferromagnetic  couplings  for  the  chains),  we  were  not  able  to  find 
this  ground  state.  Even  using  the  maximum  allowed  combinations  of  D-Wave 
function  parameters  anneal ingtime  and  numreads  (the  number  of  annealing 
runs  per  call),  the  best  that  we  could  obtain  were  solutions  corresponding  to  a 
code  size  of  7. 

On  the  other  hand,  using  a  locally  developed  embedding  solver,  which  also 
has  the  capability  to  incorporate  gauge  transformations,  we  were  able  to  find 
the  above  ground  state  corresponding  to  a  code  of  size  6,  using  the  gauge  trans¬ 
formation  shown  in  Figure  [7] 


3.3.5  Scaling  for  Larger  Cases 

While  the  B{ 2,4)  case  is  the  largest  that  can  be  solved  on  the  current  504- 
qubit  D-Wave  processor,  we  can  make  some  rough  estimates  of  how  the  qubit 
resource  requirements  scale  for  larger  cases.  For  binary  (d  =  2)  de  Bruijn 
graphs,  we  extended  the  techniques  described  above  to  generate  Ising  models 
for  the  minimum  identifying  code  problem  for  the  cases  n  =  3  through  n  =  6. 
These  models  were  then  embedded  into  ideal  Chimera  graphs  of  various  sizes 
using  the  D-Wave  software  function  sapiFindEmbeddingO .  Figure  12  shows 


how  the  number  of  qubits  needed  for  the  embedding  grows  as  a  function  of  n. 

Since  the  embedding  algorithm  is  heuristic,  it  is  possible  that  smaller  em¬ 
beddings  could  be  found,  so  the  embedding  sizes  shown  here  should  be  viewed 
as  upper  bounds.  We  computed  embeddings  firstly  for  larger  versions  of  the 
current  C8  Chimera  architecture,  which  has  8-qubit  unit  cells,  as  well  as  for 
hypothetical  C16  (16-qubit  unit  cell)  Chimera  graphs.  For  the  C8  graphs,  we 
started  with  the  current  Vesuvius  architecture,  which  contains  64  8-qubit  unit 
cells  arranged  in  an  8x8  square  grid,  and  increased  the  grid  size  to  12x12,  16x16, 
24x24,  32x32,  and  so  on  until  the  embedding  was  successful.  For  the  C16  graphs, 
we  started  with  an  8x8  grid  of  16-qubit  unit  cells,  and  increased  the  grid  size 
to  16x16,24x24,  and  32x32,  until  the  embedding  was  successful.  The  plots  show 
that  for  any  given  case,  the  embedding  size  on  C16  is  smaller  than  on  C8  due 
to  the  higher  graph  connectivity  of  the  C16  architecture;  in  other  words,  qubit 
resource  requirements  depend  on  the  hardware  graph  connectivity. 
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Figure  12:  Estimated  scaling  of  qubit  resource  requirements  for  minimum  iden¬ 
tifying  code  problem  on  binary  de  Bruijn  graphs 
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From  the  figure,  we  can  project  roughly  when  the  D-Wave  processor  would 
have  sufficient  qubits  to  accommodate  the  larger  cases  of  the  minimum  iden¬ 
tifying  code  problem  for  undirected  binary  de  Bruijn  graphs.  If  we  define  a 
processor  generation  to  be  4  times  the  number  of  qubits  as  the  previous  gener¬ 
ation  (e.g.  the  “Vesuvius”  generation  had  a  512-qubit  design  whereas  the  prior 
“Rainier”  generation  had  a  128-qubit  design),  we  can  state  that  in  roughly  1.5 
generations  we  would  be  able  to  fit  the  B( 2,  5)  case  on  the  processor.  This  is 
the  largest  d  =  2  case  for  which  the  minimum  identifying  code  size  is  presently 
known.  In  one  more  generation  beyond  that,  we  would  be  able  to  fit  the  B{ 2,6) 
case  on  the  processor,  for  which  the  minimum  identifying  code  size  is  presently 
not  known. 

Whether  these  hypothetical  D-Wave  processors  would  actually  be  able  to 
solve  these  larger  cases,  will  depend  on  the  performance  characteristics  of  those 
machines  which  is  yet  to  be  demonstrated.  The  hardware  intrinsic  control  errors 
would  need  to  be  significantly  reduced  from  the  current  levels.  For  example,  the 
embedding  we  found  for  the  £>( 2, 6)  case  on  the  C8  architecture  had  a  maximum 
chain  length  of  63  qubits.  Solving  a  problem  with  chain  lengths  this  large  would 
probably  require  6  to  7  bits  of  precision  (the  current  Vesuvius  design  has  4  bits 
of  precision).  On  the  C16  architecture,  the  embeddings  are  less  complex;  e.g. 
the  embedding  we  found  for  the  B( 2, 6)  case  only  had  a  maximum  chain  length 
of  25  qubits.  However,  it  is  not  clear  how  difficult  it  would  be  for  D-Wave  to 
achieve  a  16-qubit  unit  cell  design. 

Due  to  risk  factors  such  as  these,  it  could  take  longer  than  the  projected 
number  of  processor  generations  before  a  solution  can  be  found  to  the  B( 2,6) 
case.  On  the  other  hand,  it  may  be  possible  to  break  the  problem  down  into 
subproblems  using  the  decomposition  technique  described  above,  which  may 
make  it  easier  to  fit  the  problem  onto  the  D-Wave  processor. 

4  SMT  Solvers 

Satisfiability  Modulo  Theory  (SMT)  is  a  current  area  of  research  that  is  con¬ 
cerned  with  the  satisfiability  of  formulas  with  respect  to  some  background  the¬ 
ory  in-  SMT  solvers  combine  boolean  SAT  solving  with  decision  procedures  for 
specific  theories.  For  example,  consider  the  following  problem. 

a  =  b  +  1,  c  <  a,  c  >  b 

In  the  theory  of  the  integers,  this  problem  is  not  satisfiable  (there  are  no  integers 
a,  b,  c  where  all  the  expressions  are  true),  however  in  the  theory  of  the  real 
numbers  it  is  satisfiable  (for  example,  with  a=ll,  b=10,  c=10.5).  In  general, 
solving  an  SMT  problem  consists  of  first  solving  a  SAT  problem,  then  doing 
theory-specific  reasoning,  and  then  possibly  going  back  and  changing  the  SAT 
problem.  This  process  is  repeated  if  necessary.  In  addition,  multiple  theories  can 
also  be  used  in  the  same  SMT  problem  instance,  which  may  require  additional 
repeats  of  this  method. 
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Figure  13:  Minimum  identifying  codes  on  B(d,n) 


To  use  SMT  solvers  on  our  identifying  code  problem  for  the  undirected  de 
Bruijn  graph,  we  must  first  come  up  with  a  formulation  of  the  problem  using 
decision  procedures.  The  graph  B(d,n),  contains  dn  nodes.  For  each  of  these, 
we  create  a  boolean  variable  that  denotes  whether  or  not  the  node  is  part  of 
the  identifying  code.  We  then  also  create  an  array  of  boolean  variables  for  that 
node’s  identifying  set.  An  assertion  is  added  to  make  sure  that  each  element  of 
the  array  is  true  if  and  only  if  the  corresponding  neighbor’s  boolean  variable  is 
true  (i.e.  if  and  only  if  the  neighbor  is  part  of  the  identifying  code).  To  ensure 
unique  codes,  we  add  a  statement  to  require  that  each  node’s  identifying  set  is 
unique  from  every  other  node’s  identifying  set.  Then,  to  get  codes  of  a  fixed 
size,  we  create  an  integer  variable  for  each  node  and  add  the  constraints  that  the 
integer  is  at  least  0  and  no  greater  than  1.  Next  we  add  an  assertion  that  each 
node’s  integer  variable  is  1  if  and  only  if  its  boolean  variable  is  true.  Finally, 
we  add  a  constraint  that  the  sum  of  all  of  the  integer  variables  is  equal  to  the 
desired  identifying  code  size. 

Now  that  the  formulation  of  the  problem  has  been  determined,  we  can  use 
a  commercial  SMT  solver  to  find  solutions.  For  this  work,  we  used  the  solver 
Z3,  made  by  Microsoft  Research.  We  begin  by  first  picking  a  code  size,  and 
asking  if  there  exists  an  identifying  code  of  that  size.  If  not,  then  the  code  size 
is  increased  by  1  and  the  problem  is  posed  to  Z3  again.  This  continues  until 
an  identifying  code  of  a  specific  size  is  found.  To  find  all  satisfying  models, 
after  a  single  model  was  found  an  assertion  is  inserted  into  the  formulation 
that  requires  that  the  the  previously  found  identifying  code  be  eliminated  as 
an  option.  This  forces  Z3  to  produce  a  different  solution,  or  to  state  that  the 
formulation  is  unsatisfiable  (and  hence  no  more  identifying  codes  of  that  size 
exist).  This  process  is  repeated  in  a  loop  to  obtain  all  identifying  codes. 

Using  this  approach  on  a  single  core,  we  were  able  to  reproduce  our  results  for 
B(d,n)  from  the  HPC  method  in  much  less  time.  See  Figure  13  for  a  summary 
of  these  results.  The  numbers  in  parentheses  denote  that  we  found  a  code  of 
that  size,  but  did  not  eliminate  the  possibility  of  a  smaller  code  existing. 

Because  of  the  advancements  in  current  SAT  and  SMT  solvers,  they  offer  the 
potential  to  scale  much  better  than  a  parallelized  brute  force  approach.  This 
is  due  in  part  to  the  fact  that  many  of  today’s  solvers  are  capable  of  realizing 
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which  subsets  of  assignments  will  define  an  unsatishable  result,  and  hence  they 
will  avoid  models  in  which  those  statements  are  set.  In  our  problem,  this  might 
correspond  to  a  case  in  which  nodes  A  and  B  have  the  same  identifying  set. 
In  this  case,  the  solver  would  not  bother  looking  at  combinations  of  True/False 
assignments  on  the  other  nodes  that  do  not  affect  the  identifying  sets  of  A  or 
B. 

In  addition  to  the  sophistication  of  today’s  solvers,  there  is  also  the  possi¬ 
bility  of  parallelizing  the  search.  While  some  instances  were  run  manually  in 
a  parallel  manner  for  this  experiment,  there  is  some  research  to  be  done  on 
automatically  parallelizing  the  search  in  order  to  further  our  known  minimum 
results. 


5  Conclusions 

The  various  methods  discussed  provide  pure  mathematicians  with  a  range  of 
opportunities  for  collaboration  with  scientists  from  several  disciplines,  such  as 
computer  science  and  physics.  While  it  has  been  shown  that  some  of  these  ideas 
correspond  well  to  other  combinatorial-type  problems  (for  example,  see  0  for  a 
discussion  of  Ising  models  for  problems  from  graph  theory) ,  we  found  that  these 
approaches  alone  did  not  provide  a  strong  enough  method  for  this  problem. 
Additionally,  the  relatively  recent  appearance  of  the  D-Wave  machine  allowed 
for  very  little  literature  to  draw  from.  From  the  base  cases  that  we  constructed 
using  our  variety  of  approaches,  several  new  conjectures  have  been  developed 
and  eventually  proven  true  [2],  however  the  leg  work  needed  to  compute  the 
base  cases  required  a  deep  understanding  of  various  computing  techniques. 
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